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A theoretical study of the variation of thermodynamic and transport properties of calamitic liquid 
crystals across the isotropic-nematic phase transition is carried out by calculating the wavenumber 
(k) and time (t) dependent intermediate scattering function of the liquid, via computer simulations 
of model nematogens. The objective is to understand the experimentally observed anomalies and 
sharp variation in many thermodynamic and transport properties, namely specific heat C, sound 
attenuation coefficient Y, thermal diffusivity Dt and sound velocity c s are as the I-N transition is 
approached from the isotropic side. The small wavelength limit of the calculated intermediate scat- 
tering function F(k, t) is used to obtain the ratio of specific heats 7 and other properties mentioned 
above. We find that all of them show non-monotonic variations near the I-N transition, with F show- 
ing a cusp-like behavior. We suggest that the observed anomalous features are a direct consequence 
of the existence of pseudo-nematic domains in the system near the phase boundary and the melting 
and formation of such domains give rise to sound attenuation and also to the observed specific heat 
anomaly. A theoretical description of these anomalies should invoke translation-rotation coupling at 
molecular level. While the heterogenous dynamics observed here bear resemblance to that in deeply 
supercooled liquids near glass transition, the thermodynamic anomalies articulated here are largely 
absent in supercooled liquids. 



I. INTRODUCTION 

Time correlation functions of wavenumber dependent 
collective number density and orientational fluctuations 
in a molecular liquid play important role in understand- 
ing many microscopic aspects of dynamics in molecu- 
lar liquids. The wavenumber dependent collective cor- 
relation functions of fluctuations bridge the gap be- 
tween the single particle relaxation and the relaxation 
in the hydrodynamic limit (k — > and u> — ► O) 1 ^. 
These correlation functions can be measured in the ex- 
periments such as light scattering and inelastic neutron 
scattering£ii£i££&i2. In the hydrodynamic limit, the ex- 
pressions of these correlation functions can be obtained 
from the linearized Navier-Stokes equations 3 -^* 5 -. This 
linearized hydrodynamic model and various molecular 
extensions of them have been useful in interpreting the 
results obtained from the computer simulations and the 
experiments of simple liquids in equilibrium 3 -! 5 -^. The 
time correlation function analysis of the liquid state uses 
well-known expressions for the transport coefficients and 
the hydrodynamic behavior. Such an analysis includes 
the wavelengths and frequencies that are comparable to 
inter atomic distances and microscopic relaxation rates. 
Deviations from the hydrodynamic behavior often ap- 
pear in this analysis above some critical wavenumber. 
The wavenumbers near this critical distance may be re- 
ferred as intermediate wavenumbers. Even after many 
years of investigation by different researchers, a complete 
microscopic understanding of the relaxation of collective 
fluctuation in liquids remains incomplete. At interme- 
diate wavenumber, the correlation functions of collective 
fluctuations deviate from their hydrodynamic behavior; 
here the distance referred to is comparable to that of in- 
ter particle separation. The collective relaxation in the 
intermediate wavenumber regime can be treated in the 



framework of extended hydro dynamic s 5 ' 7 ' 1 1 ' 1 2 . 

In contrast to simple liquids, the studies of molecu- 
lar liquids are bound to be complex due to the presence 
of the internal degrees of freedom, which couple to the 
conserved modes and act as additional relaxation chan- 
nels for fluctuation. Compared to simple liquids, there 
are several parameters of intermolecular potential, like 
aspect ratio and geometry of molecule, which may affect 
the microscopic relaxation of molecular liquid s 13 ' 14 ! 15 . 

There are several attempts to study the relaxation 
of the hydrodynamic modes in molecular liquids, how- 
ever, such an analysis near the I-N transition is found 
to be absent in the literature. It is interesting to 
study this aspect of I-N transition, especially in the con- 
text of the recent experimental results by Fayer and 
coworker s 16 ! 17 ' 18 ! 19 ' 20 ! 21 ' 22 , where the orientational de- 
gree of freedom shows relaxation dynamics similar to 
that found in supercooled liquids. In another work Li 
et alw^ presented a mode coupling theory which could 
reproduce the results observed in the experiments. How- 
ever, they did not provide origin of the mode which cou- 
ples to the orientational degrees of freedom. On the 
other hand, the experimental study of collective relax- 
ation of other modes show the signature of the second 
order phase transition^. In the dynamic light scatter- 
ing experiments near the I-N transition, the spectra of 
the orientational degrees of freedom mask dynamics of 
other modes. Mode-selective dynamic light scattering 
experiments of Takagi et al^ on 5CB and 7CB liquid 
crystal samples have studied the dynamics of modes usu- 
ally hidden under fluctuation of the orientational mode. 
They have measured the thermal diffusion, orientational 
relaxation and propagating sound modes of those liq- 
uid crystal samples. Using photopyroelectric technique, 
Marinelly et. aft- have measured the thermal conduc- 
tivity and thermal diffusivity of aligned nCB samples. 
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They have found that the thermal conductivity shows 
a discontinuity near the I-N transition, while, the ther- 
mal diffusivity shows a dip. In older studies of ultrasonic 
investigations of Eden et al£2. on MBBA samples with 
out any preferred orientation with emphasis given near 
transition temperature ~ 44 °C, they found the acous- 
tic relaxation time and relaxation strength increases near 
I-N transition. 

It is now worth while to study effect of I-N transition on 
conserved degrees of freedom and also on the thermody- 
namic and transport properties. In this paper, the study 
of the collective transport properties of the nematogens 
has been restricted chiefly to the pre-transition region, 
where, even though there is an anisotropy due to the 
formation of pseudo nematic domains, they relax slowly 
to remove the anisotropy. However, as the orientational 
relaxation slows down close to the I-N transition, the re- 
laxation of the anisotropy also slows down. Since the 
I-N transition is weakly first order, there are prolonged 
pre-transitional effects, which enable the observation of 
the effect of orientational ordering in the structural re- 
laxation for a range of state points across I-N transition. 

Since I-N transition occurs due the formation of the 
orientational cage that freezes rotation of the molecules, 
as the I-N transition is approached, the exchange of the 
angular momentum with the cage gives rise to a strong 
back scattering region in the angular velocity auto corre- 
lation function (see Ref.— ). 

It is of course of interest to study what signatures 
the formation and melting of pseudo-nematic domains 
have on the transport properties, like thermal diffu- 
sivity, and also on the thermodynamic properties like 
sound velocity and specific heat. Earlier we presented 
a study of the variation of viscosity which shows weak 
anomaly across the isotropic-nematic phase transition. 
The present study is a continuation of our earlier study 26 
but now extended to include many thermodynamic and 
transport properties. Just as in the case of viscosity, it 
is not clear in the beginning what role formation of ori- 
entational order can have on properties which are largely 
determined by translational degrees of freedom. 

Investigations reported in this paper are based on the 
simulation of a model molecular liquid with aspect ratio 
~ 3 and can be approximated by prolate ellipsoids with 
two of the minor axis are equal. As reported earlie r 26 ! 27 , 
this system undergoes an I-N transition as the density 
increases from p* =0.285 to p* =0.315, gradually along 
an isotherm at T* = 1.0. The single particle relaxation 
dynamics in this system shows usual behavior as that 
of an ordinary liquid, except the presence of a marked 
anisotropy near I- N transition. 

The present study on the relaxation of fluctuation of 
collective density modes focuses on the relaxation dy- 
namics at different wavenumbers. The intermediate scat- 
tering function calculated from the simulations is fitted 
to the linearized hydrodynamic model of the intermediate 
scattering function in order to find the thermodynamic 
and transport properties of the system. Thermodynamic 



and transport properties obtained from the fitting pa- 
rameters show non-monotonic behavior at different state 
points across the I-N transition. The sound velocity ob- 
tained from the peak position of the Brillouin peak£&£ 
shows that the adiabatic sound velocity increases across 
the I-N transition; except a deviation from it's monotonic 
increase near the transition. 

An important (and we believe) new result of the 
present study is the appearance of anisotropy in the 
wavenumber dependence of the dynamic structure factor, 
F(k,t). The latter becomes a function of the direction of 
propagation of wavevector, k. 

The remaining sections of this paper are organized as 
follows. Section [II] discusses the definition of the time 
correlation functions and their expressions in the hydro- 
dynamic limit, which are used in the analysis. Details of 
simulations are presented in the section InTl The section 
IIVI presents the results and discussions on the time cor- 
relation functions calculated from the simulations. The 
summary of the work is presented in the section [Vl 



II. COLLECTIVE CORRELATION 
FUNCTIONS OF DENSITY 

For a system that consists of N linear molecules there 
are 3N coordinates of position and 2N coordinates of ori- 
entation. In this system, density is a function the posi- 
tions as well as orientations p = p(r,Sl). The number 
density of the system p(r) can be obtained from p(r, SI) 
by averaging over the orientations SI. Note that, the ori- 
entation of the molecule is not a conserved quantity. 

The intermediate structure factor F(k, SI, SI', t) which 
measures the relaxation of density fluctuations in the 
molecular liquid is defined aa 2 ^, 

F(k,n,n',t) = l( P (-k > n,o)p(k,n',i)). (i) 

This correlation function also accounts for the correla- 
tion between the orientation and density at two points 
separated in space and time. This is a generalization of 
the intermediate scattering function^' 6 -' 7 - of a monatomic 
liquid. On averaging over orientation, F(k, SI, SI' , t) gives 
F(k,t), which is the intermediate scattering function of 
number density of the system. In the theory of the lin- 
earized hydrodynamics, F(k,t) can be expressed as^ 

F(k^0,t) = S(fc^0)[ae-'/ r ° + 
(1 - a) e- t/Tb 

(cos(ujot) + b sin(u>ot))], (2) 

where S(k) is the static structure factor, a = (7 — l)/7 
(7 is the ratio of specific heats (C p /C v j), r a = 1/Dt^ 2 
(Dt is the thermal diffusivity), t/, = 1/rfc 2 (r is the 
sound wave damping constant), ujq — ck (c is the adia- 
batic sound velocity) and b = fc((3r — D v )/jc) (D v is the 
kinematic viscosity) . 
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The intermediate scattering function F(k,t) and the 
structure factor gives information about the variation of 
the thermodynamical as well as transport coefficients. 
The S(k) at the lowest value of A; gives the isothermal 
compressibility 



XT 



dp 

" ' is 



(3) 



of the system. \t of the liquid across the I-N transition 
obtained from S(k) by the relation 



Xt = — lim S(k) 



(4) 



III. DETAILS OF THE SIMULATIONS 



Molecular dynamics simulations have been carried out 
mostly for a system of 1125 Gay-Berne 29 ' 30 ellipsoids 
at temperature T*=l near I-N transition in a micro- 
Canonical ensemble. The form of the modified inter- 
molecular Gay-Berne potential used in the simulation is 

27.31.32.33.34 



U = 4e(f, Ui, Uj) 



r - a(f, Ui, Uj) + o t 
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r - a(f, Ui, Uj) + cr s 



(5) 



where f is the unit vector that passes through the center of mass of a pair of molecules, Ui and Uj unit vectors that 
passes through the major axis of a pair of ellipsoidal molecules, e and a give the strength and range of interaction, 



o(r,Ui,Uj) 



X ( (ui-r + Uj ■ f) 2 {ui 



Uj ■ ff 



1 + X{ui-Uj) l-x(ui-Uj) 



(6) 



a s is double of the minor axis b, k gives molecular elongation (aspect ratio), which is the ratio of end-to-end to 
side-to-side diameters, k = a e /a s . The aspect ratio of ellipsoids used in this simulation is 3. Here x is defined as k 
as x = (« 2 - 1)/(k 2 + 1) 



e(f, Hi, Uj) = e [1 - x 2 (ui ■ Uj)] 



V 
2 



(ili ■ f + Uj ■ f) 2 (ui ■ r - 



1 + X(ui-Uj) l-x(ui-Uj) 



(7) 



where eo is the energy parameter and x' = — 
1)/(v / k' + 1) (k' — e s /e e gives the strength of interac- 
tion which is side-to-side to end-to-end well depths). The 
value of k' used in the simulation is 5 33 i 35 i 36 , The scaling 
used for moment of inertia is /* = I /mao 2 . The density 
is scaled in the simulation as p* — pao 3 and the temper- 
ature is scaled as T* = A^T/eo- The equation of motion 
is integrated with reduced time (t* = (mcro/eo) 1 / 2 ) steps 
with At = 0.002 t*. The ellipsoid used in the simula- 
tion has minor axis b = 0.5 and major axis a = 1.5 (in 
reduced units). The simulations are done at the state 
points near the pre-transition region of phase diagram, 
shown in figure Q] The translational and rotational mo- 
tions are solved using leap-frog algorithm. The order pa- 
rameter changes dramatically in this system after density 
increases beyond 0.3. This is in accord with our previous 
works on same system—. 

The simulation starts from an equilibrated configura- 
tion of ellipsoids. Initial configuration of the ellipsoids is 
generated from a cubic lattice and then the simulation is 
run for two hundred thousand steps to obtain the equi- 
librium configuration. During the equilibration steps the 



temperature is scaled so that the system is in equilibrium 
at this particular temperature. 

In figure [1] , we show the phase diagram of ellipsoids 
with aspect ratio 3 interacting through Gay-Berne poten- 
tial. In the diagram we also show by arrow the densities 
studied at temperature (T* = 1.0). Figure [2] shows the 
variation of order parameter S — with density. 

The molecular dynamics simulations ran for 10 mil- 
lion steps to generate trajectories. We have calculated 
the collective correlation functions across I-N transition 
from this data. The wavevectors are chosen to be in the 
x, y and z directions of the simulation box. The differ- 
ence in the calculated correlation functions in different 
directions is a measure of the anisotropy developed in 
the system. Relaxation of wavenumber dependent corre- 
lation functions has been studied at three wavenumbers. 
The smallest wavenumber studied is at k — 2tt/L, where 
L is the length of the simulation box. The lowest value 
of the wavenumber k as the density increases from 0.285 
to 0.32 with a grid of 0.005 is w 0.4. The latter is studied 
with a grid of 8k w; 0.06. 
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FIG. 1: The phase diagram of the Gay-Berne ellipsoids. The 
simulations are carried out at the temperature and densities 
indicated by the arrow. 
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FIG. 2: The variation in the order parameter with density. 



IV. RESULTS AND DISCUSSIONS 
A. Density fluctuations 

Figure [3] shows the variation of the structure factor 
S(k) across the I-N transition. As density increases a 
structure of the liquid differs at different directions. In 
one of the direction (in case of present simulation) the 
director of the system is oriented approximately in the z 
direction. Note that the time scale of the director relax- 
ation is much slower than the the density relaxation- 7 . 



This enables the study of the system in different direction 
parallel and perpendicular to the director of the system. 
The development of the anisotropy in the system is not 
appreciable at the smallest wavenumbers. However, at 
large wavenumbers correlation of density increases in the 
x and y directions, at the same time along the z direc- 
tion density correlations vanish at larger wavenumbers. 
Therefore, in this system the orientational relaxation is 
mainly due to the relaxation of the local orientational 
cage rather than the global rotation of pseudo-nematic 
domain. This enables the study of the kn and k± with- 
out the use of the use of an external orienting field. Note 
that this possible only because the time scale of the were 
these correlations measured are short. 

Figure 2] shows the total intermediate scattering func- 
tion for lowest wavenumber chosen for the study of re- 
laxation, at different state points of the Gay-Berne liquid 
across I-N transition The F(k,t) at k = shows hydro- 
dynamic behavior. The rapid oscillation of the calculated 
F(k,t) at this wavenumber is related to the adiabatic 
sound propagation in the medium according to linearized 
hydrodynamics^^. As density fluctuation increases the 
short wavenumber shows a rapid growth of correlation at 
the transition point. This can be an effect of the cou- 
pling of the between the orientational and the transla- 
tional modes. This observation is not the artifact of the 
development of the global anisotropy of the system due to 
nematic domain formation. If that is the case the growth 
of correlation of the density fluctuation will be different 
in the case of x, y, and z directions of the simulation 
box. At this wavevector in all the three direction density 
relaxation shows same feature of growth of correlation 
near the transition point. Note that highest density p* 
= 0.32, which is the nematic phase, where correlation 
of static density fluctuation at the lowest wavenumber 
reduces. 



B. Thermodynamic and Transport coefficients 

F(k,t) is fitted to it's hydrodynamic form (Eq. [2] ) 
at the lowest value of k. From the parameters of the fit 
of F(k,t), at k = fci, various thermodynamic parame- 
ters and transport coefficients can be obtained. Figures 
[5] and [5] shows various thermodynamic parameters and 
transport coefficients calculated from the fit of F(k, t) for 
k = ki,. The variation of the isothermal compressibility 
of the system at I-N transition is shown in the figure^a). 
The isothermal compressibility shows a hump near I-N 
transition. The ratio of the specific heats 7 (7 = 1/1 — a) 
is plotted against density is shown in figure [5jb) . The 
variation of thermal diffusivity (Dt = {/T a k 2 ) against 
density shown in figure [HIa), which shows a dip near the 
I-N phase boundary. A similar feature was observed in 
the experiment of Marinelly et al. in the study of ther- 
mal diffusivity of aligned nCB samples^. 7 shows a rapid 
increase at the transition point. Another transport coeffi- 
cient that can be calculated from the fitting parameters is 
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FIG. 3: The calculated structure factor of the system is plotted for different densities against wavenumber. The three lines at 
each density shows the structure of the liquid directions of x, y and z: shown respectively by solid , dashed-dot and dash lines. 
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FIG. 4: The intermediate scattering function F(k, t) at lowest k plotted for different densities against time. The three lines at 
each density shows the structure of the liquid directions of x, y and z: shown respectively by solid (blue) , dashed-dot (red) 
and dash (black) lines. 



the sound wave damping constant, V (r = 1/rtfc 2 ). The 
variation of T against density is shown in figure [HJb). T 
shows a rapid growth near the TN transition. This re- 
sults of the variation sound damping constant when I-N 
transition is approached, is in qualitative agreement with 
the results of the Eden et al?^ . In figure [6fc), the adi- 
abatic sound velocity deduced from the position of the 



Brillouin peak is plotted against density. Except for a 
few points near the TN transition, the adiabatic sound 
velocity c shows a monotonic increase. In figure [6](c) , 
sound velocity shows a plateau like region very near the 
I-N phase boundary. A qualitative picture of the effects 
of the emergence of increased static correlation on the col- 
lective dynamics of a system of nematogens near the I-N 
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FIG. 5: The variation of isothermal compressibility (xt) (a) 
and the ratio of the specific heats (7) (b) are plotted against 
density across I-N transition and variation in the isothermal 
compressibility of the liquid across the I-N transition. The 
symbols diamond, circle, triangle, and left triangle represent 
respectively the data obtained from the analysis of dynamic 
intermediate scattering function along x, y, z directions and 
the average value of them. 



transition can be obtained from the analysis presented 
above. In the isotropic phase of the Gay-Berne liquid, 
the relaxation of density is isotropic due the rapid rota- 
tion of the molecules. But, near the I-N transition due to 
the freezing of individual rotation of the molecules, the 
anisotropy persists for long, which makes the density re- 
laxation anisotropic. This may be viewed in the perspec- 
tive of supercooled liquids, where it take very long time 
to relax anisotropy. In supercooled liquids there is no 
direction for the anisotropy. Near the I-N transition, in 
the isotropic liquid, the anisotropy is due to the coupling 
of number density to the transient director of the system. 
In the case of the supercooled liquids the density relax- 
ation also slows down several orders of magnitud o 37 ' 38 . 
However, in the Gay-Berne liquid the number density, 
although anisotropic, relaxes in a comparable time scale 
as that of an isotropic liquid. In the case of a liquid in 
the pre-transition region of the I-N transition the non- 
ergodicity is due the slow relaxation of the director. The 
sharp fall in 7 shows that mean square enthalpy {AH) 
fluctuations decreases compared to the same for energy 
(AE) fluctuations. If we write the enthalpy change at 
constant pressure as AH = AE + P AV , then the present 
result suggests that a positive change in AV (melting of 
pseudo-nematic domains to the isotropic phase) is accom- 
panied by a negative change in energy, which is in accor- 
dance with the understanding that the nematic phase 
has lower internal energy than the isotropic phase. On 
the other hand, the formation of a local nematic domain 
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FIG. 6: Hydrodynamic parameters obtained from the fit of 
the intermediate scattering function are plotted against den- 
sities across the I-N transition, (a) shows the variation of 
the thermal diffusivity. (b) shows the variation of the sound 
dispersion coefficient, (c) shows the variation of the sound ve- 
locity. The symbols diamond, circle, triangle, and left triangle 
represent respectively the data obtained from the analysis of 
dynamic intermediate scattering function along x, y, z direc- 
tions and the average value of them. 

due to interaction with the propagating sound wave can 
give rise to a negative AV and positive AE. Similarly, 
increase of sound attenuation is due to the interaction 
of the sound propagation with the melting/formation of 
pseudo-nematic domains which consumes energy from 
the sound wave. 



V. CONCLUDING REMARKS 

In this paper we have presented a study of I-N transi- 
tion from the perspective of molecular hydrodynamics us- 
ing molecular dynamics simulations. In a system of 1125 
particles the smallest wavenumber that can be studied 
is approximately is k ~ 0.4 in the range of the number 
densities p* = 0.285 to p* = 0.32. The calculation finds 
a fairly well defined Brillouin peak in the range of the 
wavenumbers studied. These simulations show the role 
of density fluctuation and the orientational ordering in 
development of anisotropy at intermediate wavenumebrs 
in the density relaxation near the I-N transition. How- 
ever, at large wavenumbers the calculation shows only 
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marginal evidence for the existence and propagation of 
the density fluctuations. The main results of this work 
are that it shows the effect of the formation of the director 
on the relaxation of fluctuations of the collective density 
of the system. A fit of intermediate scattering function 
to the linearized hydrodynamic model at small wavenum- 
bers (fe = 2-7T / X) show the variation of the thermody- 
namic as well as transport properties across I-N transi- 
tion. The ratio of the specific heats shows a discontinuity- 
like feature near the I-N transition. The thermal diffusion 
shows a dip; a similar behavior is reported in the recent 
experiments of Marinclly et al. 24 . The sound dispersion 
coefficient also shows a divergence like behavior near I-N 
transition. The sound velocity calculated from the peak 
position if the Brillouin peak deviate from its monotonic 
increase near I-N transition also that is in qualitative 
agreement with the experimental results of Eden^S et at. 

This work demonstrates the effect of spatial and tem- 
poral propagation of anisotropy of the molecular poten- 
tial in the medium at different state points across the I-N 
transition. Also effect of this in the thermodynamic and 
transport coefficients. Thus during the I-N transition 
the system shows the non-ergodicity due to the forma- 
tion of director, that slowly relaxes due to the mutual 
coupling of the orientation within a pseudo nematic do- 
main. Our results qualitatively agrees with experimental 
results even though the potential used is a simple ellip- 
soidal model. This shows that these results are related 
directly the anisotropy of the inter molecular potential. 

In the context of continuing interest of the comparison 
between the dynamics of supercooled liquid and liquid 
crystal o 17 ' 18 ! 27 ' 39 , it is interesting to note that the anoma- 
lies observed in sound attenuation and sound absorption 
coefficients are not known to exist (at least significantly) 
in the former. This outlines the uniqueness of the liquid 
crystalline systems where large pseudo-nematic domains 
form and melts near the I-N boundary. Such a physi- 
cal process does not exist in the supercooled liquid at a 
thermodynamic level. 

We have argued that most of the features observed can 
be qualitatively understood from the melting/formation 
of the pseudo-nematic domains. A complete theory of 



the observed hydrodynamic anomalies need to consider 
translation-rotation couplin g 40 ! 41 and is yet to be devel- 
oped. 

We are not aware of the existence of any theoretical 
treatment of the effects of slow down in the decay of 
the orientational time correlation function of nematogens 
near their I-N transition, on the isotropic dynamic struc- 
ture factor. In earlier study, we have reported observa- 
tion of an anomalous viscoelasticity of the nematogens 
near the I-N transition.- 6 " Both the problems require in- 
clusion of the effects of orientational slow down on the 
isotropic density relaxation. This is clearly a non-linear 
problem. While the orientational order parameter S is 
a non-conserved order parameter, the isotropic density p 
is a conserved order parameter. The Landau-dc Gcnncs- 
Ginzburg free energy for the isotropic phase of the liq- 
uid crystal can be written down without much difficulty. 
However, difficulty lies in the fact that the total free en- 
ergy function should be a functional of both 5(r) and 
p(r), with a coupling between the two fields. A molec- 
ular hydrodynamic theory along the above line needs to 
be developed but still appears to be non-existing. 

A valid criticism of the present work is the finite size 
of the system simulated, although N=1125 is quite large 
(and demanding) for long MD simulations, particularly 
when orientational degree of freedom is present. Also 
note that Gay-Berne is quite a complex potential to 
study. In addition to N=1125, we have also studied ex- 
tensively the size N=500. The basic features remain un- 
changed. Nevertheless, the k-^0 limit is expected to be 
sensitive to the system size, and a system size dependent 
study of the present results could be a worthwhile exer- 
cise. The intermediate wavenumber results, however, are 
not expected to be sensitive to system size. It will be 
worthwhile to study this range by neutron scattering. 
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